*******************************************************************************
********     Code to replicate Wei et al. (2024), June 12th, 2025     *********
*******************************************************************************

//This code is to present the results in our comment paper. Tables first and then figures.
// All original codes extracted from Wei et al. (2024) are marked as /*.....*/

****ado preparation
ssc install reghdfe,replace
ssc install balancetable,replace
ssc install coefplot,replace

*******************************************************************************
/* */use data.dta,clear /* */ 

/* */global controls PGDP PGDP_2 PDEN INDS CONS GOV Lnrain Lnsun /* */ 

merge 1:1 cid year using data_sup.dta,force
drop _merge

global predeter plan87 Minority special ruralpop rural2300 //selection criteria

/* */xtset cid year /* */ 

gen LnGDP=ln(GDP/deflate2018)
gen LnPGDP=ln(PGDP)
***Table 1 Diagnoses on bad controls selected by Wei et al. (2024).
foreach var in LnPGDP PDEN INDS CONS GOV Lnrain Lnsun{
	xtreg `var' NPC yr*,fe cluster(cid)
	est store m`var'
}
esttab m* , b(%9.4f) se(%7.4f) se ar2(%7.4f) nogap mtitle drop(yr*) star(* 0.1 ** 0.05 *** 0.01) 
esttab m* using Table1.rtf, replace b(%9.4f) se(%7.4f) se ar2(%7.4f) nogap drop(yr*) star(* 0.1 ** 0.05 *** 0.01) title("Table 1 Diagnoses on bad controls selected by Wei et al. (2024).")

***Table 2 Replication and examination of the core results

*column (1)
/* */xtreg CEI NPC PGDP PGDP_2 PDEN INDS CONS GOV yr*,fe cluster(cid)/* */
/* */est store m1/* */

*column (2)
/* */xtreg CEI NPC $controls yr*, fe cluster(cid)/* */
/* */est store m2/* */

*column (3)
xtreg CEI NPC yr*, fe cluster(cid)
est store m3

*column (4)
xtreg CEI NPC Lnrain Lnsun yr*, fe cluster(cid)
est store m4

*column (5)
xtreg CEI NPC PDEN INDS CONS GOV Lnrain Lnsun yr*, fe cluster(cid)
est store m5

*column (6)
xtreg CEI NPC PDEN INDS CONS Lnrain Lnsun yr*, fe cluster(cid)
est store m6

esttab m1 m2 m3 m4 m5 m6 , b(%9.4f) se(%7.4f) se ar2(%7.4f) nogap mtitle drop(yr*) star(* 0.1 ** 0.05 *** 0.01) 
esttab m1 m2 m3 m4 m5 m6 using Table2.rtf, replace b(%9.4f) se(%7.4f) se ar2(%7.4f) nogap drop(yr*) star(* 0.1 ** 0.05 *** 0.01) title("Table 2 Replication and examination of the core results")

***Table 3 Replication and examination of the core results

/* */forvalu x=2010/2018{
       preserve
             keep if year==`x'
                 set seed 0001
             psmatch2 treat $controls, outcome(CEI) n(5) ate ties caliper(0.05) common quietly
             drop if _weight==.
             save  psmdid`x'.dta,replace
       restore
}
//here please make a pause and rerun the following append code to get the matched sample
 preserve 
use psmdid2010.dta ,clear     
forvalu x=2011/2018{
    append using psmdid`x'.dta   
}/* */

*column (1)
xtreg CEI NPC $controls yr*,fe cluster(cid) 
est store psm1

*column (2)
reghdfe CEI NPC $controls [aw = _weight],a(cid year) cluster(cid) keepsingletons
est store psm2 

*column (3)
reghdfe CEI NPC [aw = _weight],a(cid year) cluster(cid) keepsingletons
est store psm3   

esttab psm1 psm2 psm3 , b(%9.4f) se(%7.4f) se ar2(%7.4f) nogap mtitle drop(yr*) star(* 0.1 ** 0.05 *** 0.01) 
   
restore //here please make a pause and rerun this code to return back the basic data

***Table 4 Balance test.

*Panel A
preserve 
keep if year ==2010
balancetable treat $controls using "Table 4.xlsx", vce(robust) format(%8.4f) noobservations starlevels(* 0.1 ** 0.05 *** 0.01) ctitles("Treatment group" "Control group" "Unconditional diff.") sheet(Sheet1) modify

*Panel B
balancetable treat $predeter using "Table 4.xlsx", vce(robust) format(%8.4f) noobservations starlevels(* 0.1 ** 0.05 *** 0.01) ctitles("Treatment group" "Control group" "Unconditional diff.") sheet(Sheet2) modify

restore

***Table 5 Estimated results of augmented DID specification.

*column (1)
xtreg CEI NPC c.post2014#c.($predeter) yr*, fe cluster(cid)
est store c1 

*column (2)
gen t_1 = year-2010+1
gen t_2 = t_1^2
gen t_3 = t_1^3
xtreg CEI NPC c.(t_1 t_2 t_3)#c.($predeter) yr*, fe cluster(cid)
est store c2 

*column (3)
xtreg CEI NPC i.year#c.($predeter) yr*, fe cluster(cid)
est store c3 
esttab c1 c2 c3 , b(%9.4f) se(%7.4f) se ar2(%7.4f) nogap mtitle keep(NPC) star(* 0.1 ** 0.05 *** 0.01)

***Table 6 Estimated results of augmented DID specification.

*column (1)
xtreg LnGDP NPC i.year#c.($predeter) yr*, fe cluster(cid)
est store g1 

*column (2)
/* */gen event=year-2014 if treat==1
tab event,gen(eventz)
forvalue i=1/9{
          replace eventz`i'=0 if eventz`i'==.

}
drop eventz1/* */
xtreg LnGDP eventz* i.year#c.($predeter) yr*, fe cluster(cid)
test eventz2 = eventz3 = eventz4 = 0
est store g2 
esttab g1 g2, b(%9.4f) se(%7.4f) se ar2(%7.4f) nogap mtitle keep(NPC eventz*) star(* 0.1 ** 0.05 *** 0.01)


***Figure 2. Event study estimates.

*Panel A
/* */xtreg CEI eventz* $controls yr*,fe cluster(cid)  

coefplot , coeflabels (eventz2 = "-3"  eventz3 = "-2"  eventz4 = "-1"  eventz5 = "0"   ///
    eventz6 = "1"   ///
    eventz7 = "2"    ///
    eventz8 = "3"    ///
    eventz9 = "4" )    ///
    vertical baselevels keep(*eventz*) yline(0) levels(95)  addplot(line @b @at) ///
	    ciopts(lpattern(dash) recast(rcap) msize(medium)) ///
		 msymbol(circle_hollow)   ///
       scheme(s1mono)  ///
	   yline(0, lwidth(vthin) lpattern(dash) lcolor(blace)) ///
	   xtitle("time", size(small))  ///
	   ylabel(-0.15(0.05)0.15,format(%5.2f) )/* */  

*Panel B
xtreg CEI eventz* yr*,fe cluster(cid)  

coefplot , coeflabels (eventz2 = "-3"  eventz3 = "-2"  eventz4 = "-1"  eventz5 = "0"   ///
    eventz6 = "1"   ///
    eventz7 = "2"    ///
    eventz8 = "3"    ///
    eventz9 = "4" )    ///
    vertical baselevels keep(*eventz*) yline(0) levels(95)  addplot(line @b @at) ///
	    ciopts(lpattern(dash) recast(rcap) msize(medium)) ///
		 msymbol(circle_hollow)   ///
       scheme(s1mono)  ///
	   yline(0, lwidth(vthin) lpattern(dash) lcolor(blace)) ///
	   xtitle("time", size(small))  ///
	   ylabel(-0.15(0.05)0.15,format(%5.2f) ) 	   

*Panel C
xtreg CEI eventz* PDEN INDS CONS Lnrain Lnsun yr*,fe cluster(cid)  

coefplot , coeflabels (eventz2 = "-3"  eventz3 = "-2"  eventz4 = "-1"  eventz5 = "0"   ///
    eventz6 = "1"   ///
    eventz7 = "2"    ///
    eventz8 = "3"    ///
    eventz9 = "4" )    ///
    vertical baselevels keep(*eventz*) yline(0) levels(95)  addplot(line @b @at) ///
	    ciopts(lpattern(dash) recast(rcap) msize(medium)) ///
		 msymbol(circle_hollow)   ///
       scheme(s1mono)  ///
	   yline(0, lwidth(vthin) lpattern(dash) lcolor(blace)) ///
	   xtitle("time", size(small))  ///
	   ylabel(-0.15(0.05)0.15,format(%5.2f) ) 	   

*Panel D
xtreg CEI eventz* i.year#c.($predeter) yr*,fe cluster(cid)  

coefplot , coeflabels (eventz2 = "-3"  eventz3 = "-2"  eventz4 = "-1"  eventz5 = "0"   ///
    eventz6 = "1"   ///
    eventz7 = "2"    ///
    eventz8 = "3"    ///
    eventz9 = "4" )    ///
    vertical baselevels keep(*eventz*) yline(0) levels(95)  addplot(line @b @at) ///
	    ciopts(lpattern(dash) recast(rcap) msize(medium)) ///
		 msymbol(circle_hollow)   ///
       scheme(s1mono)  ///
	   yline(0, lwidth(vthin) lpattern(dash) lcolor(blace)) ///
	   xtitle("time", size(small))  ///
	   ylabel(-0.15(0.05)0.15,format(%5.2f) ) 	   
	      

	   
	   